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Exotic hadronic states and all-to-all quark propagators 

UKQCD Collaboration: Chris Michael and Janne Peisa a * 

a Theoretical Physics Division, Department of Mathematical Sciences, 
University of Liverpool, Liverpool L69 3BX, United Kingdom 

We discuss methods to obtain accurate hadronic spectra with propagating quarks. Comparing the determination 
of masses for spin-exotic hybrid mesons with glueball mass determinations, we conclude that quark propagators 
from all sites to all other sites would enable great improvement in the errors. Such propagators are achievable 
by using stochastic estimators. We discuss previous attempts and present our method for maximal variance 
reduction. This is a very promising technique and we illustrate it by obtaining the spectrum of ground state and 
excited B mesons in the limit where the b quark is static. 



On 



> 

m 

o 
in 
o 
t> 

On 
-t— > 

Jh 

Oh! 



X 



1. INTRODUCTION 

One of the goals of lattice QCD is to pin down 
the hadronic spectrum accurately in the quenched 
approximation and then to explore the shifts as 
the dynamical quark degrees of freedom are re- 
introduced. This is particularly relevant to guide 
experimental searches for glueballs and hybrid 
mesons. The situation for lattice determination 
of the spectrum of hybrid mesons, at present, is 
that the first clear signals have been obtained but 
much greater precision would be valuable. 

It is interesting to contrast this situation with 
the glueball mass determination for which a pre- 
cise continuum value has been obtained in 
the quenched approximation. The glueball cor- 
relator is measured from all sites on the lattice 
which gives large statistics even from one gauge 
configuration. For correlations involving gauge 
links only, this is very easy to achieve computa- 
tionally. For correlations involving fermion prop- 
agation, the inversion of the fermion matrix from 
even a single source is computationally expen- 
sive. Thus the hybrid correlator has been mea- 
sured ||,||] from only a single source for each 
gauge configuration which gives a rather noisy 
signal. Moreover to construct appropriate non- 
local sources for hybrid meson studies can also 
involve extra quark propagator inversions which 
is costly. The same constraints also apply to a 
detailed study of the orbital and radial excited 
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Figure 1. The effective mass of the J PC — 1 + 
hybrid meson versus t from the UKQCD data |] 
with U-shaped source and several different sinks. 



hadrons as well. 

The result of the recent UKQCD determina- 
tion H of the spin-exotic J PC = 1 h hybrid 
meson is shown in Figure 1. This was obtained 
from 350 propagators from each of two sources 
separated by a fixed distance and joined by a U- 
shaped gauge link. By using different operators 
at the sink, several different correlators were eval- 
uated which allows a separation of ground state 
and excited states of the required quantum num- 
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bers. The study with a clover- Wilson fermion ac- 
tion using Csw — 1-4785 at = 6.0 and a hop- 
ping parameter corresponding to a quark mass 
close to strange, gave a hybrid meson mass of 
2.0(0.2) GeV for the ss state. The results for 
other exotic spin-parity hybrid mesons are con- 
sistent with being heavier than this 1 _+ state. 

As already indicated, a more precise determi- 
nation of these spectra would be valuable. With 
this in mind, we now discuss the feasibility of eval- 
uating all-to-all hadronic correlators. As an illus- 
tration of the potency of this approach, we will 
present a preliminary study of the spectrum of ex- 
cited bq mesons. Here spin-exotic hybrid mesons 
cannot be identified because charge conjugation 
C is not a good quantum number. Nevertheless, 
lattice studies can guide experimental searches in 
this area where few experimental results are avail- 
able. 

2. STOCHASTIC PROPAGATORS 

The conventional approach to quark propaga- 
tion in a gauge background field (this discussion 
applies equally to the quenched case and to full 
QCD with dynamical quarks) is to solve the lat- 
tice Dirac equation Qip = 5 iteratively from a 
single source to determine a propagator G sat- 
isfying QG = 1. This problem is equivalent to 
that of the inversion of a sparse matrix since Q 
has only local and nearest neighbour contribu- 
tions - this inversion is solved by iteration which 
is computationally expensive. For a local source 
at x, the propagator G(y, x) contains information 
about the quark propagation from this source to 
all lattice sites y. For the Wilson-Dirac discreti- 
sation (this includes the SW-clover case) then 
Q(x,y) = 75Q(y, a;)*75 and hence the propaga- 
tion from all sites to the source is also obtained. 
Thus in the usual approach, only propagators to 
or from x = (0, 0, 0, 0) are determined. This lim- 
its the correlations that can be constructed to ex- 
plore hadronic spectra and matrix elements and 
can cause large statistical errors and large correla- 
tions among different t values since so few samples 
of the correlation are evaluated for each gauge 
configuration. 

In the case of a meson with one static quark 



and one propagating light quark, then only the 
elements of G with the same spatial component 
x = y will be used in such correlations and two 
such numbers (from positive and negative time 
directions) per gauge configuration will be avail- 
able for use to study the correlation versus time 
t. Thus only a very small part of the information 
obtained labouriously by inverting Q to obtain G 
is used in this extreme case. An attractive al- 
ternative is to obtain the propagator from many 
different sites. The default way to achieve this is 
just to repeat the inversion for sources at differ- 
ent positions. This is computationally demanding 
and, if the propagators were to be stored, would 
involve prodigious storage requirements for all-to- 
all propagators. 

An alternative route is to accept stochastic esti- 
mators of the propagators. Since the fermion ma- 
trix inversion is conducted in a gauge field sam- 
ple, it is not necessary to determine the propa- 
gator G to 10 digit accuracy as is usually done. 
Indeed, varying the elements of the propagator 
randomly by a few percent makes little change to 
the hadronic correlators obtained. The problem is 
to achieve this stochastic sampling in an unbiased 
way. Just stopping the inversion algorithm at rel- 
ative accuracy of 1% introduces bias. One accept- 
able approach is to use a pseudo-fermion method 
with scalar fields interacting via the fermion ma- 
trix. This amounts to obtaining a stochastic es- 
timate of the inverse of a matrix M by use of the 
gaussian Monte Carlo process: 

Z = j Dcj> e -^' M * (1) 

with M^ 1 —< 4>*4>j >■ This approach requires 
a positive definite matrix M which can be ar- 
ranged by taking M = Q^Q where Q is the 
Wilson-Dirac fermion matrix. Then M will be 
sparse (it has at most next to nearest neighbour 
terms) and the Monte Carlo process can be sim- 
ulated efficiently. Taking a = 1 , ... N samples 
of the scalar fields <f>f yields stochastic estimates 
of the inverse of M. Thence the propagator 
Q^ 1 can be recovered from these samples since 

Gji = Qji = < {Qik<t>k)*<f>j >■ 

The stochastic evaluation of propagators has 
the advantage that the number of correlations 
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measured increases hugely - by L 3 T for a lattice of 
spatial size L and time length T. This is compen- 
sated by the increase in error from the stochastic 
(rather than exact) propagators used. The sta- 
tistical error is easily estimated since each scalar 
field <pi has a variance of order 1 (here we nor- 
malise the Wilson-Dirac matrix as Q = 1 — KD 
with hopping parameter K) . Thus for N samples 
each propagator will have errors of order TV -1 / 2 . 
This will result in a statistical error on hadron 
correlations which is independent of time separa- 
tion t. So for small t one will have small relative 
errors while at larger t the relative error will be 
huge since the signal is so small. Since hadron 
correlators at reasonably large t values are needed 
to separate out ground states and excited states, 
a way to reduce these statistical errors is essential. 

One suggestion || is to use a local multi-hit 
variance reduction for the <p fields. This is equiv- 
alent to taking an improved estimator for the 
<p field which is the average over many Monte 
Carlo updates of that field with its neighbours 
held fixed. This is very easy to implement and 
provides a significant error reduction. This ap- 
proach has been used for B-meson physics suc- 
cessfully g. Such a variance reduction is a wel- 
come amelioration but still does not give an ab- 
solute error decreasing with t. This very attrac- 
tive goal can, however, be achieved and here we 
present our proposals. 

Imagine that instead of a multi-hit update of 
one site in the presence of its neighbours held 
fixed, one considered the updating of a fixed re- 
gion P. Let the scalar fields within P be labelled 
<p and those on the boundary of P be labelled s. 
Then the required improved estimate of <p with s 
held fixed is given by 

Vi = jp J D(p(p i e-^^ lCl ^ k+ ^^ kSk+h ^ (2) 

where we have distinguished the elements of M 
connecting the set of (p fields inside P to them- 
selves (M) and those connecting them to the 
boundary (M). The integral over <p is gaussian 
and we obtain: 

Vi = -Mf/MjkSk (3) 
This can be visualised as the propagator M _1 in 



region P from source Ms at the boundary to site 
i. We will call v the maximally variance reduced 
stochastic estimator. If site i lies a minimum of 
d links from the nearest boundary point, then we 
would expect u, to decrease in magnitude as d 
increases. In the lowest order hopping parameter 
expansion this behaviour will be as K . Since the 
variance of s will remain of order one, this implies 
that the variance of v will be maximally reduced 
by taking a partition P which has a boundary 
as far as possible from the point i. However, we 
must remember that this variance reduction for- 
mula applies if just one <p field within P is to have 
its variance reduced. To evaluate a variance re- 
duced propagator one must choose two disjoint 
regions P and R and solve for the variance re- 
duced fields Vi and Wj in P and R respectively. 
Then we will have Gji = < v*Wj > with no bias. 
Thus the rules of operation are to take regions P 
and R around points i and j such that the bound- 
aries of P and R are as far from these points as 
possible. 

One geometry which allows this is to take P as 
the time slices < t < T/2 and R as the time 
slices T/2 < t < T. The boundary sources S are 
then the time slices at t = and T/2. Then prop- 
agators between P and R (also between P and S 
and between R and S) can be evaluated with max- 
imal variance reduction. This is not all-to-all just 
most-to-most, although by using in turn several 
different choices of partitioning into P and R, all- 
to-all can be obtained. We have achieved the goal 
of stochastic propagators from most sites to most 
sites but at what computational cost? Each vari- 
ance reduction operation corresponds to a sparse 
matrix inversion with a source Ms (where s are 
the stochastic samples of (p on the source planes 
at t = and T/2) in a region less than half the 
lattice size. N=12 samples of v and w can be 
obtained for a computational effort equivalent to 
one usual inversion of M for all 12 colour-spin 
components at one source. This increases by a 
factor of 2 or so the computation of evaluating 
by Monte Carlo the stochastic fields in the first 
place. This is more than balanced by the huge 
reduction in error as we shall see. 

Since this procedure is conceptually a triply 
nested Monte Carlo, we summarise the steps 
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needed. 

• Create gauge configurations g. 

• In create sample scalar fields for all x and 
colour-spin components by Monte Carlo us- 
ing equation 1: <j) a (g). 

• Using the scalar field as a source on time- 
planes at t = and T/2, solve iteratively 
using equation 3 for the improved estima- 
tors v a (g) and w a (g). 

. then Gji = E 9 H a (Q^{g)Yw«(g). 

This procedure will certainly provide maxi- 
mally variance reduced stochastic estimators of 
propagators. The number of the stochastic sam- 
ples N of the scalar field should be chosen so that 
the resulting error on the quantity of interest from 
one gauge configuration is small enough but not 
so small that it is less than the intrinsic variation 
from one gauge configuration to another. We ex- 
plore first the case of the B meson using a static 
approximation for the heavy quark. This is a very 
favourable case for stochastic inversion since there 
is a large increase in statistics (through using all 
space points as sources) and the hadronic correla- 
tion of interest is proportional to the light quark 
propagator so the choice of N is not critical. Our 
choice of partitions P and R is also seen to be ap- 
propriate because propagators are only needed in 
the time direction for which the variance reduc- 
tion will have maximal effect. 

3. B SPECTROSCOPY 

We consider a small lattice (8 3 16) at = 5.74 
with Wilson fermions of hopping parameter K = 
0.156. This choice was motivated by pre-existing 
studies ||. As a first example we evaluated the 
B meson correlator at time separation t using lo- 
cal hadronic operators at source and sink. Then 
we compared conventional inversion with various 
implementations of stochastic inversion. Sample 
results are shown in Table 1 where the compari- 
son has been made for equal disk storage of prop- 
agators or scalars. We chose N = 25 samples to 
evaluate the stochastic estimates of the propaga- 
tors. In the maximally reduced variance method, 
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Figure 2. The B meson effective mass versus 
t from our data with different combinations of 
local, fuzzed and jiDi sources and sinks together 
with a three exponential fit. Also shown (squares) 
are the Wuppertal data || for smeared source and 
local sink from 170 propagators. 



Vi and Wj were chosen symmetrically either side 
of the source time planes. 

The maximal variance reduction gives a fac- 
tor of 7 improvement in error compared to other 
methods for only an overall computational in- 
crease of a factor of 4. This is equivalent to 
a net gain of a factor of 12 in computing time 
for a similar result. Moreover, the stochas- 
tic method allows correlations involving different 
sources (smeared, fuzzed, etc) to be constructed 
at little extra cost. This can be seen in Fig- 
ure 2 which shows the effective mass plots for 
some of these correlators - these results come 
from stochastic evaluation from 20480 sources - 
as described in Table 1. The results [|| from con- 
ventional inversions (with 170 sources) are seen to 
be significantly less precise than those obtained 
here. 

Also illustrated is a three state fit to our 
hadronic correlations which yields ma = 0.855(6) 
for the ground state and m'a — 1.233(15) for 
the first excited state of the same quantum num- 
bers, where the errors quoted are statistical only. 
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Table 1 

B meson correlators at t = 7. 



Method 


C(7) x 10 Y 


Data Set 


CPU 


MR inversion 


3712(147) 


propagators from 4 sources 


1 






for 10 gauge fields 




Stochastic inversion 




25 samples of </> 




Basic 


2754(926) 


for 20 gauge fields 


2 


Local multihit 


3418(410) 




2 


Maximal variance reduction 


3761(21) 
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We have also extracted the correlations corre- 
sponding to the L-excited states (P and D- 
wave). Of particular interest is the splitting 
of the P wave level into two energy values in 
the static heavy quark limit, corresponding to 
J p = + , 2 + mesons. Our preliminary results 
are ma — 1.35(5) and 1.24(4) respectively, where 
the errors are purely statistical. Note that, if the 
level ordering indicated by these preliminary re- 
sults should be substantiated, this is a first indi- 
cation of a meson system for which the 2 + state 
is lighter than the + . This is not unexpected 
in potential based models since the relevant spin- 
orbit splitting has this sign at large distance |^,[| . 
Further work using larger lattices, several lattice 
spacings and several hopping parameters is under 
way to explore this further. 

4. HADRONIC OBSERVABLES 

For observables with more than one light quark 
propagator, it is possible to use the stochastic 
method, but some care is needed. One way to 
grasp the subtlety is to imagine that there are 
two quarks with different flavours and to split 
the sum over samples a = 1 , . . . N into subsets 
for each flavour. Provided these two subsets are 
independent, then propagators of each flavour are 
obtained without bias. In practice, provided the 
stochastic estimators 4> a are independent of a, 
then a sum over a, (3 = 1, . . .N with a ^ (3 
gives the required information. This illustrates 
that the number of samples is effectively N 2 for 
mesons and N 3 for baryons. Indeed we have ob- 
served that the error from the stochastic sampling 
reduces as 1/N for correlators involving two light 
quarks in the region of small N. For large N val- 



ues, there will be no advantage in reducing the 
error from one gauge sample to below the intrin- 
sic variation between gauge samples. 

We discuss briefly the potential advantages 
of using stochastic propagators for different 
hadronic observables. As we have shown above, 
for heavy quarks in the static limit, light quark 
propagators are needed in the time direction and 
these are efficiently produced by the stochastic 
method with maximal variance reduction. Thus 
other static observables such as Ah and the matrix 
element Bb are efficiently evaluated. 

For mesons made of light quarks, it is neces- 
sary to project the correlations into a specific mo- 
mentum state. This implies a sum of correlations 
over the whole spatial volume at that time. Even 
though the signal will be localised at small dif- 
ference between spatial coordinates of source and 
sink, in this sum the noise will come from all spa- 
tial sites. Another way to see this is that the 
propagator for t time steps and x space steps will 
behave in the lowest order hopping parameter ex- 
pansion as K t+X whereas the variance reduction 
we have employed above will reduce the error by 
K l only. Thus it appears that stochastic inversion 
will require a relatively large number of samples 
N to reduce this noise in the case of the momen- 
tum zero hadronic correlation. 

One interesting area of application is to discon- 
nected quark diagrams - as needed for rj meson 
studies and glueball mixing. The maximal vari- 
ance reduction technique does not seem promising 
since the variance reduced average of 4> 2 will in- 
volve Mj" 1 w hich is n °t easily obtained computa- 
tionally (this is after all the origin of the problem 
here anyway). It is feasible to use a more local 
region P and this has been studied R]. 
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5. DISCUSSION 

Particularly for dynamical fermion configura- 
tions which are computationally expensive to ob- 
tain, there is a virtue in measuring the hadronic 
properties from a gauge configuration as com- 
pletely as possible. We have presented a method 
to extract correlators from all sources to all sinks 
to achieve this. Moreover the storage require- 
ments are not excessive. For the case of correla- 
tors involving one static quark, the gain in signal 
is very significant. For exploration of observables 
involving more than one source and one sink, such 
as matrix elements or meson-meson interactions, 
then our proposal is very promising. 
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